Shells

Init symbols for sympy


In [1]:
from sympy import *
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
x1, x2, x3 = symbols("x_1 x_2 x_3", real = True)
alpha1, alpha2, alpha3 = symbols("alpha_1 alpha_2 alpha_3", real = True, positive=True)
R, L, ga, gv = symbols("R L g_a g_v", real = True, positive=True)
init_printing()

Cylindrical coordinates


In [2]:
a1 = pi / 2 + (L / 2 - alpha1)/R

x = R * cos(a1)
y = alpha2
z = R * sin(a1)

r = x*N.i + y*N.j + z*N.k

Mid-surface coordinates is defined with the following vector $\vec{r}=\vec{r}(\alpha_1, \alpha_2)$


In [3]:
r


Out[3]:
$$(- R \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\alpha_{2})\mathbf{\hat{j}_{N}} + (R \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

Tangent to curve


In [4]:
r1 = r.diff(alpha1)
r2 = r.diff(alpha2)
k1 = trigsimp(r1.magnitude())
k2 = trigsimp(r2.magnitude())
r1 = r1/k1 
r2 = r2/k2

In [5]:
r1


Out[5]:
$$(\cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

In [6]:
r2


Out[6]:
$$\mathbf{\hat{j}_{N}}$$

Normal to curve


In [7]:
n = r1.cross(r2)
n = trigsimp(n.normalize())
n


Out[7]:
$$(- \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

Base Vectors $\vec{R}_1, \vec{R}_2, \vec{R}_3$


In [8]:
R_alpha=r+alpha3*n
R_alpha


Out[8]:
$$(- R \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - \alpha_{3} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\alpha_{2})\mathbf{\hat{j}_{N}} + (R \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + \alpha_{3} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

In [9]:
R1=R_alpha.diff(alpha1)
R2=R_alpha.diff(alpha2)
R3=R_alpha.diff(alpha3)

In [10]:
trigsimp(R1)


Out[10]:
$$(\left(1 + \frac{\alpha_{3}}{R}\right) \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\left(1 + \frac{\alpha_{3}}{R}\right) \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

In [11]:
R2


Out[11]:
$$\mathbf{\hat{j}_{N}}$$

In [12]:
R3


Out[12]:
$$(- \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

Draw


In [13]:
%matplotlib inline

import matplotlib.pyplot as plt
import plot

x = (R + alpha3) * cos(a1)
z = (R + alpha3) * sin(a1)

alpha1_x = lambdify([R, L, alpha1, alpha3], x, "numpy")
alpha3_z = lambdify([R, L, alpha1, alpha3], z, "numpy")

R_num = 1/0.8
L_num = 2

x1_start = 0
x1_end = L_num
x3_start = -0.05
x3_end = 0.05

def alpha_to_x(a1, a2, a3):
    x=alpha1_x(R_num, L_num, a1, a3)
    z=alpha3_z(R_num, L_num, a1, a3)
    return x, 0, z
    
plot.plot_init_geometry_2(x1_start, x1_end, x3_start, x3_end, alpha_to_x)


Base Vectors $\vec{R}^1, \vec{R}^2, \vec{R}^3$


In [14]:
eps=trigsimp(R1.dot(R2.cross(R3)))
R_1=simplify(trigsimp(R2.cross(R3)/eps))
R_2=simplify(trigsimp(R3.cross(R1)/eps))
R_3=simplify(trigsimp(R1.cross(R2)/eps))
eps


Out[14]:
$$1 + \frac{\alpha_{3}}{R}$$

In [15]:
R_1


Out[15]:
$$(\frac{R}{R + \alpha_{3}} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\frac{R}{R + \alpha_{3}} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

In [16]:
R_2


Out[16]:
$$\mathbf{\hat{j}_{N}}$$

In [17]:
R_3


Out[17]:
$$(- \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{i}_{N}} + (\cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )})\mathbf{\hat{k}_{N}}$$

Jacobi matrix:

$ A = \left( \begin{array}{ccc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} & \frac{\partial x_1}{\partial \alpha_3} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \frac{\partial x_3}{\partial \alpha_1} & \frac{\partial x_3}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \end{array} \right)$

$ \left[ \begin{array}{ccc} \vec{R}_1 & \vec{R}_2 & \vec{R}_3 \end{array} \right] = \left[ \begin{array}{ccc} \vec{e}_1 & \vec{e}_2 & \vec{e}_3 \end{array} \right] \cdot \left( \begin{array}{ccc} \frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} & \frac{\partial x_1}{\partial \alpha_3} \\ \frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \frac{\partial x_3}{\partial \alpha_1} & \frac{\partial x_3}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \\ \end{array} \right) = \left[ \begin{array}{ccc} \vec{e}_1 & \vec{e}_2 & \vec{e}_3 \end{array} \right] \cdot A$

$ \left[ \begin{array}{ccc} \vec{e}_1 & \vec{e}_2 & \vec{e}_3 \end{array} \right] =\left[ \begin{array}{ccc} \vec{R}_1 & \vec{R}_2 & \vec{R}_3 \end{array} \right] \cdot A^{-1}$


In [18]:
dx1da1=R1.dot(N.i)
dx1da2=R2.dot(N.i)
dx1da3=R3.dot(N.i)

dx2da1=R1.dot(N.j)
dx2da2=R2.dot(N.j)
dx2da3=R3.dot(N.j)

dx3da1=R1.dot(N.k)
dx3da2=R2.dot(N.k)
dx3da3=R3.dot(N.k)

A=Matrix([[dx1da1, dx1da2, dx1da3], [dx2da1, dx2da2, dx2da3], [dx3da1, dx3da2, dx3da3]])
simplify(A)


Out[18]:
$$\left[\begin{matrix}\frac{1}{R} \left(R + \alpha_{3}\right) \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & - \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\\0 & 1 & 0\\\frac{1}{R} \left(R + \alpha_{3}\right) \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\end{matrix}\right]$$

In [19]:
A_inv = trigsimp(A**-1)
simplify(trigsimp(A_inv))


Out[19]:
$$\left[\begin{matrix}\frac{R}{R + \alpha_{3}} \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & \frac{R}{R + \alpha_{3}} \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\\0 & 1 & 0\\- \sin{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & \cos{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\end{matrix}\right]$$

In [20]:
trigsimp(A.det())


Out[20]:
$$1 + \frac{\alpha_{3}}{R}$$

Metric tensor

${\displaystyle \hat{G}=\sum_{i,j} g^{ij}\vec{R}_i\vec{R}_j}$


In [21]:
g11=R_1.dot(R_1)
g12=R_1.dot(R_2)
g13=R_1.dot(R_3)

g21=R_2.dot(R_1)
g22=R_2.dot(R_2)
g23=R_2.dot(R_3)

g31=R_3.dot(R_1)
g32=R_3.dot(R_2)
g33=R_3.dot(R_3)

G=Matrix([[g11, g12, g13],[g21, g22, g23], [g31, g32, g33]])
G=trigsimp(G)
G


Out[21]:
$$\left[\begin{matrix}\frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$

${\displaystyle \hat{G}=\sum_{i,j} g_{ij}\vec{R}^i\vec{R}^j}$


In [22]:
g_11=R1.dot(R1)
g_12=R1.dot(R2)
g_13=R1.dot(R3)

g_21=R2.dot(R1)
g_22=R2.dot(R2)
g_23=R2.dot(R3)

g_31=R3.dot(R1)
g_32=R3.dot(R2)
g_33=R3.dot(R3)

G_con=Matrix([[g_11, g_12, g_13],[g_21, g_22, g_23], [g_31, g_32, g_33]])
G_con=trigsimp(G_con)
G_con


Out[22]:
$$\left[\begin{matrix}\frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$

In [23]:
G_inv = G**-1
G_inv


Out[23]:
$$\left[\begin{matrix}\frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$$

Christoffel symbols


In [24]:
DIM = 3

G_con_diff = MutableDenseNDimArray.zeros(DIM, DIM, DIM)
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            xdiff = alpha1
            if (k == 0): 
                xdiff = alpha1
            elif (k == 1):
                xdiff = alpha2
            elif (k == 2):
                xdiff = alpha3
            
            G_con_diff[i,j,k]=G_con[i,j].diff(xdiff)
            


GK = MutableDenseNDimArray.zeros(DIM, DIM, DIM)
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            res = S(0)
            for m in range(DIM):
                res = res + G[m,k]*(G_con_diff[i,m,j]+G_con_diff[j,m,i]-G_con_diff[i,j,m])
            GK[i,j,k] = simplify(S(1)/S(2)*res)


GK


Out[24]:
$$\left[\begin{matrix}\left[\begin{matrix}0 & 0 & - \frac{1}{R^{2}} \left(R + \alpha_{3}\right)\\0 & 0 & 0\\\frac{1}{R + \alpha_{3}} & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}\frac{1}{R + \alpha_{3}} & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\end{matrix}\right]$$

Gradient of vector

$ \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_3 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \\ \nabla_3 u_2 \\ \nabla_1 u_3 \\ \nabla_2 u_3 \\ \nabla_3 u_3 \\ \end{array}

\right)

B \cdot \left( \begin{array}{c} u_1 \\ \frac { \partial u_1 } { \partial \alpha_1} \\ \frac { \partial u_1 } { \partial \alpha_2} \\ \frac { \partial u_1 } { \partial \alpha_3} \\ u_2 \\ \frac { \partial u_2 } { \partial \alpha_1} \\ \frac { \partial u_2 } { \partial \alpha_2} \\ \frac { \partial u_2 } { \partial \alpha_3} \\ u_3 \\ \frac { \partial u_3 } { \partial \alpha_1} \\ \frac { \partial u_3 } { \partial \alpha_2} \\ \frac { \partial u_3 } { \partial \alpha_3} \\ \end{array} \right) = B \cdot D \cdot \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ \frac { \partial u^1 } { \partial \alpha_3} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \frac { \partial u^2 } { \partial \alpha_3} \\ u^3 \\ \frac { \partial u^3 } { \partial \alpha_1} \\ \frac { \partial u^3 } { \partial \alpha_2} \\ \frac { \partial u^3 } { \partial \alpha_3} \\ \end{array} \right) $


In [60]:
def row_index_to_i_j_grad(i_row):
    return i_row // 3, i_row % 3
        

B = zeros(9, 12)
B[0,1] = S(1)
B[1,2] = S(1)

B[2,3] = S(1)

B[3,5] = S(1)
B[4,6] = S(1)
B[5,7] = S(1)

B[6,9] = S(1)
B[7,10] = S(1)
B[8,11] = S(1)

for row_index in range(9):
    i,j=row_index_to_i_j_grad(row_index)
    B[row_index, 0] = -GK[i,j,0]
    B[row_index, 4] = -GK[i,j,1]
    B[row_index, 8] = -GK[i,j,2]

B


Out[60]:
$$\left[\begin{array}{cccccccccccc}0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{1}{R + \alpha_{3}} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\- \frac{1}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

In [26]:
q=(1+alpha3/R)
g_=(q*q)
koef=g_.diff(alpha3)

In [27]:
D=eye(12)
D[0,0]=G_con[0,0]
D[1,1]=G_con[0,0]
D[2,2]=G_con[0,0]
D[3,3]=G_con[0,0]
D[3,0]=G_con[0,0].diff(alpha3)
D


Out[27]:
$$\left[\begin{array}{cccccccccccc}\frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\frac{1}{R^{2}} \left(2 R + 2 \alpha_{3}\right) & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

In [28]:
simplify(B*D)


Out[28]:
$$\left[\begin{array}{cccccccccccc}0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & 0\\0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\frac{1}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\- \frac{1}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

Strain tensor

$ \left( \begin{array}{c} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{array}

\right)

E \cdot \left( \begin{array}{c} \nabla_1 u_1 \\ \nabla_2 u_1 \\ \nabla_3 u_1 \\ \nabla_1 u_2 \\ \nabla_2 u_2 \\ \nabla_3 u_2 \\ \nabla_1 u_3 \\ \nabla_2 u_3 \\ \nabla_3 u_3 \\ \end{array} \right)$


In [29]:
E=zeros(6,9)
E[0,0]=1
E[1,4]=1
E[2,8]=1
E[3,1]=1
E[3,3]=1
E[4,2]=1
E[4,6]=1
E[5,5]=1
E[5,7]=1
E


Out[29]:
$$\left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0\end{matrix}\right]$$

In [30]:
Q=E*B*D
Q=simplify(Q)

Q


Out[30]:
$$\left[\begin{array}{cccccccccccc}0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0\end{array}\right]$$

In [31]:
def E_NonLinear(grad_u):
    N = 3

    du = zeros(N, N)

    #    print("===Deformations===")

    for i in range(N):
        for j in range(N):
            index = i*N+j
            du[j,i] = grad_u[index]

    #    print("========")

    a_values = S(1)/S(2) * du * G


    E_NL = zeros(6,9)
    E_NL[0,0] = a_values[0,0]
    E_NL[0,3] = a_values[0,1]
    E_NL[0,6] = a_values[0,2]

    E_NL[1,1] = a_values[1,0]
    E_NL[1,4] = a_values[1,1]
    E_NL[1,7] = a_values[1,2]

    E_NL[2,2] = a_values[2,0]
    E_NL[2,5] = a_values[2,1]
    E_NL[2,8] = a_values[2,2]

    E_NL[3,1] = 2*a_values[0,0]
    E_NL[3,4] = 2*a_values[0,1]
    E_NL[3,7] = 2*a_values[0,2]

    E_NL[4,0] = 2*a_values[2,0]
    E_NL[4,3] = 2*a_values[2,1]
    E_NL[4,6] = 2*a_values[2,2]

    E_NL[5,2] = 2*a_values[1,0]
    E_NL[5,5] = 2*a_values[1,1]
    E_NL[5,8] = 2*a_values[1,2]


    return E_NL

u1 = Function("u_1")
u2 = Function("u_2")
u3 = Function("u_3")

gu = zeros(12,1) 
gu[0] = u1(alpha1,alpha2,alpha3)
gu[1] = u1(alpha1,alpha2,alpha3).diff(alpha1)
gu[2] = u1(alpha1,alpha2,alpha3).diff(alpha2)
gu[3] = u1(alpha1,alpha2,alpha3).diff(alpha3)

gu[4] = u2(alpha1,alpha2,alpha3)
gu[5] = u2(alpha1,alpha2,alpha3).diff(alpha1)
gu[6] = u2(alpha1,alpha2,alpha3).diff(alpha2)
gu[7] = u2(alpha1,alpha2,alpha3).diff(alpha3)

gu[8] = u3(alpha1,alpha2,alpha3)
gu[9] = u3(alpha1,alpha2,alpha3).diff(alpha1)
gu[10] = u3(alpha1,alpha2,alpha3).diff(alpha2)
gu[11] = u3(alpha1,alpha2,alpha3).diff(alpha3)

gradu=B*D*gu


E_NL = E_NonLinear(gradu)
E_NL


Out[31]:
$$\left[\begin{matrix}\frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(\frac{\left(R + \alpha_{3}\right)^{2}}{2 R^{2}} \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} + \frac{\operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}}{2 R^{2}} \left(R + \alpha_{3}\right)\right) & 0 & 0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{2}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} - \frac{\operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}}{2 R^{2}} \left(R + \alpha_{3}\right) & 0 & 0\\0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{2}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{2}} \operatorname{u_{2}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{2}} \operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0\\0 & 0 & \frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(\frac{1}{2} \left(- \frac{1}{R^{2}} \left(R + \alpha_{3}\right) + \frac{1}{R^{2}} \left(2 R + 2 \alpha_{3}\right)\right) \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} + \frac{\left(R + \alpha_{3}\right)^{2}}{2 R^{2}} \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}\right) & 0 & 0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{2}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{1}{2} \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}\\0 & \frac{2 R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(\frac{\left(R + \alpha_{3}\right)^{2}}{2 R^{2}} \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} + \frac{\operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}}{2 R^{2}} \left(R + \alpha_{3}\right)\right) & 0 & 0 & \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{2}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} - \frac{1}{R^{2}} \left(R + \alpha_{3}\right) \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0\\\frac{2 R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(\frac{1}{2} \left(- \frac{1}{R^{2}} \left(R + \alpha_{3}\right) + \frac{1}{R^{2}} \left(2 R + 2 \alpha_{3}\right)\right) \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} + \frac{\left(R + \alpha_{3}\right)^{2}}{2 R^{2}} \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}\right) & 0 & 0 & \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{2}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0\\0 & 0 & \frac{\partial}{\partial \alpha_{2}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{\partial}{\partial \alpha_{2}} \operatorname{u_{2}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )} & 0 & 0 & \frac{\partial}{\partial \alpha_{2}} \operatorname{u_{3}}{\left (\alpha_{1},\alpha_{2},\alpha_{3} \right )}\end{matrix}\right]$$

Elasticity tensor(stiffness tensor)

General form


In [32]:
from sympy import MutableDenseNDimArray
C_x = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
                el = Symbol(elem_index)
                C_x[i,j,k,l] = el
                
C_x


Out[32]:
$$\left[\begin{matrix}\left[\begin{matrix}C^{1111} & C^{1112} & C^{1113}\\C^{1121} & C^{1122} & C^{1123}\\C^{1131} & C^{1132} & C^{1133}\end{matrix}\right] & \left[\begin{matrix}C^{1211} & C^{1212} & C^{1213}\\C^{1221} & C^{1222} & C^{1223}\\C^{1231} & C^{1232} & C^{1233}\end{matrix}\right] & \left[\begin{matrix}C^{1311} & C^{1312} & C^{1313}\\C^{1321} & C^{1322} & C^{1323}\\C^{1331} & C^{1332} & C^{1333}\end{matrix}\right]\\\left[\begin{matrix}C^{2111} & C^{2112} & C^{2113}\\C^{2121} & C^{2122} & C^{2123}\\C^{2131} & C^{2132} & C^{2133}\end{matrix}\right] & \left[\begin{matrix}C^{2211} & C^{2212} & C^{2213}\\C^{2221} & C^{2222} & C^{2223}\\C^{2231} & C^{2232} & C^{2233}\end{matrix}\right] & \left[\begin{matrix}C^{2311} & C^{2312} & C^{2313}\\C^{2321} & C^{2322} & C^{2323}\\C^{2331} & C^{2332} & C^{2333}\end{matrix}\right]\\\left[\begin{matrix}C^{3111} & C^{3112} & C^{3113}\\C^{3121} & C^{3122} & C^{3123}\\C^{3131} & C^{3132} & C^{3133}\end{matrix}\right] & \left[\begin{matrix}C^{3211} & C^{3212} & C^{3213}\\C^{3221} & C^{3222} & C^{3223}\\C^{3231} & C^{3232} & C^{3233}\end{matrix}\right] & \left[\begin{matrix}C^{3311} & C^{3312} & C^{3313}\\C^{3321} & C^{3322} & C^{3323}\\C^{3331} & C^{3332} & C^{3333}\end{matrix}\right]\end{matrix}\right]$$

Include symmetry


In [33]:
C_x_symmetry = MutableDenseNDimArray.zeros(3, 3, 3, 3)

def getCIndecies(index):
    if (index == 0):
        return 0, 0
    elif (index == 1):
        return 1, 1
    elif (index == 2):
        return 2, 2
    elif (index == 3):
        return 0, 1
    elif (index == 4):
        return 0, 2
    elif (index == 5):
        return 1, 2
    
for s in range(6):
    for t in range(s, 6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
        el = Symbol(elem_index)
        C_x_symmetry[i,j,k,l] = el
        C_x_symmetry[i,j,l,k] = el
        C_x_symmetry[j,i,k,l] = el
        C_x_symmetry[j,i,l,k] = el
        C_x_symmetry[k,l,i,j] = el
        C_x_symmetry[k,l,j,i] = el
        C_x_symmetry[l,k,i,j] = el
        C_x_symmetry[l,k,j,i] = el

                
C_x_symmetry


Out[33]:
$$\left[\begin{matrix}\left[\begin{matrix}C^{1111} & C^{1112} & C^{1113}\\C^{1112} & C^{1122} & C^{1123}\\C^{1113} & C^{1123} & C^{1133}\end{matrix}\right] & \left[\begin{matrix}C^{1112} & C^{1212} & C^{1213}\\C^{1212} & C^{2212} & C^{1223}\\C^{1213} & C^{1223} & C^{3312}\end{matrix}\right] & \left[\begin{matrix}C^{1113} & C^{1213} & C^{1313}\\C^{1213} & C^{2213} & C^{1323}\\C^{1313} & C^{1323} & C^{3313}\end{matrix}\right]\\\left[\begin{matrix}C^{1112} & C^{1212} & C^{1213}\\C^{1212} & C^{2212} & C^{1223}\\C^{1213} & C^{1223} & C^{3312}\end{matrix}\right] & \left[\begin{matrix}C^{1122} & C^{2212} & C^{2213}\\C^{2212} & C^{2222} & C^{2223}\\C^{2213} & C^{2223} & C^{2233}\end{matrix}\right] & \left[\begin{matrix}C^{1123} & C^{1223} & C^{1323}\\C^{1223} & C^{2223} & C^{2323}\\C^{1323} & C^{2323} & C^{3323}\end{matrix}\right]\\\left[\begin{matrix}C^{1113} & C^{1213} & C^{1313}\\C^{1213} & C^{2213} & C^{1323}\\C^{1313} & C^{1323} & C^{3313}\end{matrix}\right] & \left[\begin{matrix}C^{1123} & C^{1223} & C^{1323}\\C^{1223} & C^{2223} & C^{2323}\\C^{1323} & C^{2323} & C^{3323}\end{matrix}\right] & \left[\begin{matrix}C^{1133} & C^{3312} & C^{3313}\\C^{3312} & C^{2233} & C^{3323}\\C^{3313} & C^{3323} & C^{3333}\end{matrix}\right]\end{matrix}\right]$$

Isotropic material


In [34]:
C_isotropic = MutableDenseNDimArray.zeros(3, 3, 3, 3)

C_isotropic_matrix = zeros(6)

mu = Symbol('mu')
la = Symbol('lambda')

for s in range(6):
    for t in range(s, 6):
        if (s < 3 and t < 3):
            if(t != s):
                C_isotropic_matrix[s,t] = la
                C_isotropic_matrix[t,s] = la
            else:
                C_isotropic_matrix[s,t] = 2*mu+la
                C_isotropic_matrix[t,s] = 2*mu+la
        elif (s == t):
            C_isotropic_matrix[s,t] = mu
            C_isotropic_matrix[t,s] = mu
            
for s in range(6):
    for t in range(s, 6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        el = C_isotropic_matrix[s, t]
        C_isotropic[i,j,k,l] = el
        C_isotropic[i,j,l,k] = el
        C_isotropic[j,i,k,l] = el
        C_isotropic[j,i,l,k] = el
        C_isotropic[k,l,i,j] = el
        C_isotropic[k,l,j,i] = el
        C_isotropic[l,k,i,j] = el
        C_isotropic[l,k,j,i] = el

                
C_isotropic


Out[34]:
$$\left[\begin{matrix}\left[\begin{matrix}\lambda + 2 \mu & 0 & 0\\0 & \lambda & 0\\0 & 0 & \lambda\end{matrix}\right] & \left[\begin{matrix}0 & \mu & 0\\\mu & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & \mu\\0 & 0 & 0\\\mu & 0 & 0\end{matrix}\right]\\\left[\begin{matrix}0 & \mu & 0\\\mu & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}\lambda & 0 & 0\\0 & \lambda + 2 \mu & 0\\0 & 0 & \lambda\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & \mu\\0 & \mu & 0\end{matrix}\right]\\\left[\begin{matrix}0 & 0 & \mu\\0 & 0 & 0\\\mu & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & \mu\\0 & \mu & 0\end{matrix}\right] & \left[\begin{matrix}\lambda & 0 & 0\\0 & \lambda & 0\\0 & 0 & \lambda + 2 \mu\end{matrix}\right]\end{matrix}\right]$$

In [35]:
def getCalpha(C, A, q, p, s, t):
    res = S(0)
    for i in range(3):
        for j in range(3):        
            for k in range(3):
                for l in range(3):
                    res += C[i,j,k,l]*A[q,i]*A[p,j]*A[s,k]*A[t,l]
    return simplify(trigsimp(res))
                    


C_isotropic_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                c = getCalpha(C_isotropic, A_inv, i, j, k, l)
                C_isotropic_alpha[i,j,k,l] = c

C_isotropic_alpha[0,0,0,0]


Out[35]:
$$\frac{R^{4} \left(\lambda + 2 \mu\right)}{\left(R + \alpha_{3}\right)^{4}}$$

In [36]:
C_isotropic_matrix_alpha = zeros(6)

for s in range(6):
    for t in range(6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        C_isotropic_matrix_alpha[s,t] = C_isotropic_alpha[i,j,k,l]
        
C_isotropic_matrix_alpha


Out[36]:
$$\left[\begin{matrix}\frac{R^{4} \left(\lambda + 2 \mu\right)}{\left(R + \alpha_{3}\right)^{4}} & \frac{R^{2} \lambda}{\left(R + \alpha_{3}\right)^{2}} & \frac{R^{2} \lambda}{\left(R + \alpha_{3}\right)^{2}} & 0 & 0 & 0\\\frac{R^{2} \lambda}{\left(R + \alpha_{3}\right)^{2}} & \lambda + 2 \mu & \lambda & 0 & 0 & 0\\\frac{R^{2} \lambda}{\left(R + \alpha_{3}\right)^{2}} & \lambda & \lambda + 2 \mu & 0 & 0 & 0\\0 & 0 & 0 & \frac{R^{2} \mu}{\left(R + \alpha_{3}\right)^{2}} & 0 & 0\\0 & 0 & 0 & 0 & \frac{R^{2} \mu}{\left(R + \alpha_{3}\right)^{2}} & 0\\0 & 0 & 0 & 0 & 0 & \mu\end{matrix}\right]$$

Orthotropic material


In [37]:
C_orthotropic = MutableDenseNDimArray.zeros(3, 3, 3, 3)

C_orthotropic_matrix = zeros(6)

for s in range(6):
    for t in range(s, 6):
        elem_index = 'C^{{{}{}}}'.format(s+1, t+1)
        el = Symbol(elem_index)
        if ((s < 3 and t < 3) or t == s):
            C_orthotropic_matrix[s,t] = el
            C_orthotropic_matrix[t,s] = el
            
for s in range(6):
    for t in range(s, 6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        el = C_orthotropic_matrix[s, t]
        C_orthotropic[i,j,k,l] = el
        C_orthotropic[i,j,l,k] = el
        C_orthotropic[j,i,k,l] = el
        C_orthotropic[j,i,l,k] = el
        C_orthotropic[k,l,i,j] = el
        C_orthotropic[k,l,j,i] = el
        C_orthotropic[l,k,i,j] = el
        C_orthotropic[l,k,j,i] = el

                
C_orthotropic


Out[37]:
$$\left[\begin{matrix}\left[\begin{matrix}C^{11} & 0 & 0\\0 & C^{12} & 0\\0 & 0 & C^{13}\end{matrix}\right] & \left[\begin{matrix}0 & C^{44} & 0\\C^{44} & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & C^{55}\\0 & 0 & 0\\C^{55} & 0 & 0\end{matrix}\right]\\\left[\begin{matrix}0 & C^{44} & 0\\C^{44} & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}C^{12} & 0 & 0\\0 & C^{22} & 0\\0 & 0 & C^{23}\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & C^{66}\\0 & C^{66} & 0\end{matrix}\right]\\\left[\begin{matrix}0 & 0 & C^{55}\\0 & 0 & 0\\C^{55} & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & C^{66}\\0 & C^{66} & 0\end{matrix}\right] & \left[\begin{matrix}C^{13} & 0 & 0\\0 & C^{23} & 0\\0 & 0 & C^{33}\end{matrix}\right]\end{matrix}\right]$$

Orthotropic material in shell coordinates


In [38]:
def getCalpha(C, A, q, p, s, t):
    res = S(0)
    for i in range(3):
        for j in range(3):        
            for k in range(3):
                for l in range(3):
                    res += C[i,j,k,l]*A[q,i]*A[p,j]*A[s,k]*A[t,l]
    return simplify(trigsimp(res))
                    


C_orthotropic_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                c = getCalpha(C_orthotropic, A_inv, i, j, k, l)
                C_orthotropic_alpha[i,j,k,l] = c

C_orthotropic_alpha[0,0,0,0]


Out[38]:
$$\frac{R^{4}}{\left(R + \alpha_{3}\right)^{4}} \left(C^{11} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right)$$

In [39]:
C_orthotropic_matrix_alpha = zeros(6)

for s in range(6):
    for t in range(6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        C_orthotropic_matrix_alpha[s,t] = C_orthotropic_alpha[i,j,k,l]
        
C_orthotropic_matrix_alpha


Out[39]:
$$\left[\begin{matrix}\frac{R^{4}}{\left(R + \alpha_{3}\right)^{4}} \left(C^{11} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & \frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(C^{12} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & \frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(C^{11} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & 0 & \frac{R^{3}}{8 \left(R + \alpha_{3}\right)^{3}} \left(- 2 C^{11} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - C^{11} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 2 C^{13} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 2 C^{33} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - C^{33} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 4 C^{55} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )}\right) & 0\\\frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(C^{12} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & C^{22} & C^{12} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & - \frac{R \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )}}{2 R + 2 \alpha_{3}} \left(C^{12} - C^{23}\right) & 0\\\frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(C^{11} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & C^{12} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & C^{11} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & \frac{R}{8 \left(R + \alpha_{3}\right)} \left(- 2 C^{11} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + C^{11} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - 2 C^{13} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 2 C^{33} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + C^{33} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - 4 C^{55} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )}\right) & 0\\0 & 0 & 0 & \frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(C^{44} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{66} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & 0 & - \frac{R \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )}}{2 R + 2 \alpha_{3}} \left(C^{44} - C^{66}\right)\\\frac{R^{3}}{8 \left(R + \alpha_{3}\right)^{3}} \left(- 2 C^{11} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - C^{11} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 2 C^{13} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 2 C^{33} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - C^{33} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 4 C^{55} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )}\right) & - \frac{R \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )}}{2 R + 2 \alpha_{3}} \left(C^{12} - C^{23}\right) & \frac{R}{8 \left(R + \alpha_{3}\right)} \left(- 2 C^{11} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + C^{11} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - 2 C^{13} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + 2 C^{33} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + C^{33} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - 4 C^{55} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )}\right) & 0 & \frac{R^{2}}{\left(R + \alpha_{3}\right)^{2}} \left(C^{11} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{55} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 2 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{55} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\right) & 0\\0 & 0 & 0 & - \frac{R \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )}}{2 R + 2 \alpha_{3}} \left(C^{44} - C^{66}\right) & 0 & C^{44} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{66} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\end{matrix}\right]$$

Physical coordinates

$u^1=\frac{u_{[1]}}{1+\frac{\alpha_3}{R}}$

$\frac{\partial u^1} {\partial \alpha_3}=\frac{1}{1+\frac{\alpha_3}{R}} \frac{\partial u_{[1]}} {\partial \alpha_3} + u_{[1]} \frac{\partial} {\partial \alpha_3} \left( \frac{1}{1+\frac{\alpha_3}{R}} \right) =\frac{1}{1+\frac{\alpha_3}{R}} \frac{\partial u_{[1]}} {\partial \alpha_3} - u_{[1]} \frac{1}{R \left( 1+\frac{\alpha_3}{R} \right)^2} $


In [40]:
H1=simplify(sqrt(G_con[0,0]))
H2=simplify(sqrt(G_con[1,1]))
H3=simplify(sqrt(G_con[2,2]))

P=zeros(12,12)
P[0,0]=1/H1
P[1,0]=(1/H1).diff(alpha1)
P[1,1]=1/H1
P[2,0]=(1/H1).diff(alpha2)
P[2,2]=1/H1
P[3,0]=(1/H1).diff(alpha3)
P[3,3]=1/H1

P[4,4]=1/H2
P[5,4]=(1/H2).diff(alpha1)
P[5,5]=1/H2
P[6,4]=(1/H2).diff(alpha2)
P[6,6]=1/H2
P[7,4]=(1/H2).diff(alpha3)
P[7,7]=1/H2

P[8,8]=1/H3
P[9,8]=(1/H3).diff(alpha1)
P[9,9]=1/H3
P[10,8]=(1/H3).diff(alpha2)
P[10,10]=1/H3
P[11,8]=(1/H3).diff(alpha3)
P[11,11]=1/H3
P=simplify(P)
P


Out[40]:
$$\left[\begin{array}{cccccccccccc}\frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{R}{\left(R + \alpha_{3}\right)^{2}} & 0 & 0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

In [41]:
Grad_U = simplify(B*D*P)
rows, cols = Grad_U.shape
Grad_U_P=zeros(rows, cols)

for i in range(3):
    for j in range(3):
        ratio=1
        if (i==0):
            ratio = ratio*H1
        elif (i==1):
            ratio = ratio*H2
        elif (i==2):
            ratio = ratio*H3
            
        if (j==0):
            ratio = ratio*H1
        elif (j==1):
            ratio = ratio*H2
        elif (j==2):
            ratio = ratio*H3
            
        
        row_index = i*3+j
        for column_index in range(cols):
            Grad_U_P[row_index, column_index]=Grad_U[row_index, column_index]/ratio

Grad_U_P = simplify(Grad_U_P)
Grad_U_P


Out[41]:
$$\left[\begin{array}{cccccccccccc}0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{R + \alpha_{3}} & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\- \frac{1}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{R}{R + \alpha_{3}} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

In [42]:
StrainL=simplify(E*Grad_U_P)
StrainL


Out[42]:
$$\left[\begin{array}{cccccccccccc}0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{R + \alpha_{3}} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & 0 & 1 & 0 & 0 & \frac{R}{R + \alpha_{3}} & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{1}{R + \alpha_{3}} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & \frac{R}{R + \alpha_{3}} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0\end{array}\right]$$

In [43]:
u1 = Function("u_1")
u2 = Function("u_2")
u3 = Function("u_3")

gu = zeros(12,1) 
gu[0] = u1(alpha1,alpha3)
gu[1] = u1(alpha1,alpha3).diff(alpha1)
gu[2] = u1(alpha1,alpha3).diff(alpha2)
gu[3] = u1(alpha1,alpha3).diff(alpha3)


gu[8] = u3(alpha1)
gu[9] = u3(alpha1).diff(alpha1)


gradup=Grad_U_P*gu


StrainNL = E_NonLinear(gradup)*gradup
simplify(StrainNL)


Out[43]:
$$\left[\begin{matrix}\frac{1}{2 \left(R + \alpha_{3}\right)^{4}} \left(R^{2} \left(R \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{3} \right )} + \operatorname{u_{3}}{\left (\alpha_{1} \right )}\right)^{2} + \left(R + \alpha_{3}\right)^{2} \left(R \frac{d}{d \alpha_{1}} \operatorname{u_{3}}{\left (\alpha_{1} \right )} - \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{3} \right )}\right)^{2}\right)\\0\\\frac{R^{2} \left(\frac{\partial}{\partial \alpha_{3}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{3} \right )}\right)^{2}}{2 \left(R + \alpha_{3}\right)^{2}}\\0\\\frac{R^{2}}{\left(R + \alpha_{3}\right)^{3}} \left(R \frac{\partial}{\partial \alpha_{1}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{3} \right )} + \operatorname{u_{3}}{\left (\alpha_{1} \right )}\right) \frac{\partial}{\partial \alpha_{3}} \operatorname{u_{1}}{\left (\alpha_{1},\alpha_{3} \right )}\\0\end{matrix}\right]$$

Stiffness tensor


In [44]:
C_isotropic_alpha_p = MutableDenseNDimArray.zeros(3, 3, 3, 3)
q=1+alpha3/R
for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                fact = 1
                if (i==0):
                    fact = fact*q
                if (j==0):
                    fact = fact*q
                if (k==0):
                    fact = fact*q
                if (l==0):
                    fact = fact*q
                C_isotropic_alpha_p[i,j,k,l] = simplify(C_isotropic_alpha[i,j,k,l]*fact)
            
C_isotropic_matrix_alpha_p = zeros(6)

for s in range(6):
    for t in range(6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        C_isotropic_matrix_alpha_p[s,t] = C_isotropic_alpha_p[i,j,k,l]
        
C_isotropic_matrix_alpha_p


Out[44]:
$$\left[\begin{matrix}\lambda + 2 \mu & \lambda & \lambda & 0 & 0 & 0\\\lambda & \lambda + 2 \mu & \lambda & 0 & 0 & 0\\\lambda & \lambda & \lambda + 2 \mu & 0 & 0 & 0\\0 & 0 & 0 & \mu & 0 & 0\\0 & 0 & 0 & 0 & \mu & 0\\0 & 0 & 0 & 0 & 0 & \mu\end{matrix}\right]$$

In [45]:
C_orthotropic_alpha_p = MutableDenseNDimArray.zeros(3, 3, 3, 3)
q=1+alpha3/R
for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                fact = 1
                if (i==0):
                    fact = fact*q
                if (j==0):
                    fact = fact*q
                if (k==0):
                    fact = fact*q
                if (l==0):
                    fact = fact*q
                C_orthotropic_alpha_p[i,j,k,l] = simplify(C_orthotropic_alpha[i,j,k,l]*fact)
            
C_orthotropic_matrix_alpha_p = zeros(6)

for s in range(6):
    for t in range(6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        C_orthotropic_matrix_alpha_p[s,t] = C_orthotropic_alpha_p[i,j,k,l]
        
C_orthotropic_matrix_alpha_p


Out[45]:
$$\left[\begin{matrix}C^{11} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & C^{12} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & C^{11} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & - \frac{C^{11}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{11}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{13}}{4} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{33}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{33}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{55}}{2} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} & 0\\C^{12} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & C^{22} & C^{12} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & \frac{1}{2} \left(- C^{12} + C^{23}\right) \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} & 0\\C^{11} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{13} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & C^{12} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{23} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & C^{11} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + 4 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & - \frac{C^{11}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{11}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{13}}{4} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{33}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{33}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{55}}{2} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} & 0\\0 & 0 & 0 & C^{44} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{66} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0 & \frac{1}{2} \left(- C^{44} + C^{66}\right) \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )}\\- \frac{C^{11}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{11}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{13}}{4} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{33}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{33}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{55}}{2} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} & \frac{1}{2} \left(- C^{12} + C^{23}\right) \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} & - \frac{C^{11}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{11}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{13}}{4} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{33}}{4} \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} + \frac{C^{33}}{8} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} - \frac{C^{55}}{2} \sin{\left (\frac{2}{R} \left(L - 2 \alpha_{1}\right) \right )} & 0 & C^{11} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 2 C^{13} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{33} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{55} \sin^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} - 2 C^{55} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{55} \cos^{4}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} & 0\\0 & 0 & 0 & \frac{1}{2} \left(- C^{44} + C^{66}\right) \sin{\left (\frac{1}{R} \left(L - 2 \alpha_{1}\right) \right )} & 0 & C^{44} \sin^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )} + C^{66} \cos^{2}{\left (\frac{1}{R} \left(\frac{L}{2} - \alpha_{1}\right) \right )}\end{matrix}\right]$$

Square of segment

$A=\frac {\theta}{2} \left( R + h_2 \right)^2-\frac {\theta}{2} \left( R + h_1 \right)^2$


In [46]:
theta, h1, h2=symbols('theta h_1 h_2')
square_geom=theta/2*(R+h2)**2-theta/2*(R+h1)**2
expand(simplify(square_geom))


Out[46]:
$$- R h_{1} \theta + R h_{2} \theta - \frac{h_{1}^{2} \theta}{2} + \frac{h_{2}^{2} \theta}{2}$$

${\displaystyle A=\int_{0}^{L}\int_{h_1}^{h_2} \left( 1+\frac{\alpha_3}{R} \right) d \alpha_1 d \alpha_3}, L=R \theta$


In [47]:
square_int=integrate(integrate(1+alpha3/R, (alpha3, h1, h2)), (alpha1, 0, theta*R))
expand(simplify(square_int))


Out[47]:
$$- R h_{1} \theta + R h_{2} \theta - \frac{h_{1}^{2} \theta}{2} + \frac{h_{2}^{2} \theta}{2}$$

Tymoshenko theory

$u^1 \left( \alpha_1, \alpha_2, \alpha_3 \right)=u\left( \alpha_1 \right)+\alpha_3\gamma \left( \alpha_1 \right) $

$u^2 \left( \alpha_1, \alpha_2, \alpha_3 \right)=0 $

$u^3 \left( \alpha_1, \alpha_2, \alpha_3 \right)=w\left( \alpha_1 \right) $

$ \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ \frac { \partial u^1 } { \partial \alpha_3} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \frac { \partial u^2 } { \partial \alpha_3} \\ u^3 \\ \frac { \partial u^3 } { \partial \alpha_1} \\ \frac { \partial u^3 } { \partial \alpha_2} \\ \frac { \partial u^3 } { \partial \alpha_3} \\ \end{array} \right) = T \cdot \left( \begin{array}{c} u \\ \frac { \partial u } { \partial \alpha_1} \\ \gamma \\ \frac { \partial \gamma } { \partial \alpha_1} \\ w \\ \frac { \partial w } { \partial \alpha_1} \\ \end{array} \right) $


In [48]:
T=zeros(12,6)
T[0,0]=1
T[0,2]=alpha3
T[1,1]=1
T[1,3]=alpha3
T[3,2]=1

T[8,4]=1
T[9,5]=1
T


Out[48]:
$$\left[\begin{matrix}1 & 0 & \alpha_{3} & 0 & 0 & 0\\0 & 1 & 0 & \alpha_{3} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 1\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$

In [49]:
D_p_T = StrainL*T
K = Symbol('K')
D_p_T = D_p_T.subs(R, 1/K)
simplify(D_p_T)


Out[49]:
$$\left[\begin{matrix}0 & \frac{1}{K \alpha_{3} + 1} & 0 & \frac{\alpha_{3}}{K \alpha_{3} + 1} & \frac{K}{K \alpha_{3} + 1} & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & \frac{1}{K \alpha_{3} + 1} & 0 & 0 & \frac{1}{K \alpha_{3} + 1}\\0 & 0 & 0 & 0 & 0 & 0\end{matrix}\right]$$

In [50]:
u = Function("u")
t = Function("theta")
w = Function("w")

u1=u(alpha1)+alpha3*t(alpha1)
u3=w(alpha1)

gu = zeros(12,1) 
gu[0] = u1
gu[1] = u1.diff(alpha1)
gu[3] = u1.diff(alpha3)


gu[8] = u3
gu[9] = u3.diff(alpha1)


gradup=Grad_U_P*gu

simplify(gradup.subs(R,1/K))

# o20=(K*u(alpha1)-w(alpha1).diff(alpha1)+t(alpha1))/2
# o21=K*t(alpha1)
# O=1/2*o20*o20+alpha3*o20*o21-alpha3*K/2*o20*o20
# O=expand(O)
# O=collect(O,alpha3)
# simplify(O)

StrainNL = E_NonLinear(gradup)*gradup
StrainNL = StrainNL.subs(R,1/K)
simplify(StrainNL)


Out[50]:
$$\left[\begin{matrix}\frac{1}{2 \left(K \alpha_{3} + 1\right)^{4}} \left(\left(K \alpha_{3} + 1\right)^{2} \left(K \left(\alpha_{3} \theta{\left (\alpha_{1} \right )} + u{\left (\alpha_{1} \right )}\right) - \frac{d}{d \alpha_{1}} w{\left (\alpha_{1} \right )}\right)^{2} + \left(K w{\left (\alpha_{1} \right )} + \alpha_{3} \frac{d}{d \alpha_{1}} \theta{\left (\alpha_{1} \right )} + \frac{d}{d \alpha_{1}} u{\left (\alpha_{1} \right )}\right)^{2}\right)\\0\\\frac{\theta^{2}{\left (\alpha_{1} \right )}}{2 \left(K \alpha_{3} + 1\right)^{2}}\\0\\\frac{\theta{\left (\alpha_{1} \right )}}{\left(K \alpha_{3} + 1\right)^{3}} \left(K w{\left (\alpha_{1} \right )} + \alpha_{3} \frac{d}{d \alpha_{1}} \theta{\left (\alpha_{1} \right )} + \frac{d}{d \alpha_{1}} u{\left (\alpha_{1} \right )}\right)\\0\end{matrix}\right]$$

Square theory

$u^1 \left( \alpha_1, \alpha_2, \alpha_3 \right)=u_{10}\left( \alpha_1 \right)p_0\left( \alpha_3 \right)+u_{11}\left( \alpha_1 \right)p_1\left( \alpha_3 \right)+u_{12}\left( \alpha_1 \right)p_2\left( \alpha_3 \right) $

$u^2 \left( \alpha_1, \alpha_2, \alpha_3 \right)=0 $

$u^3 \left( \alpha_1, \alpha_2, \alpha_3 \right)=u_{30}\left( \alpha_1 \right)p_0\left( \alpha_3 \right)+u_{31}\left( \alpha_1 \right)p_1\left( \alpha_3 \right)+u_{32}\left( \alpha_1 \right)p_2\left( \alpha_3 \right) $

$ \left( \begin{array}{c} u^1 \\ \frac { \partial u^1 } { \partial \alpha_1} \\ \frac { \partial u^1 } { \partial \alpha_2} \\ \frac { \partial u^1 } { \partial \alpha_3} \\ u^2 \\ \frac { \partial u^2 } { \partial \alpha_1} \\ \frac { \partial u^2 } { \partial \alpha_2} \\ \frac { \partial u^2 } { \partial \alpha_3} \\ u^3 \\ \frac { \partial u^3 } { \partial \alpha_1} \\ \frac { \partial u^3 } { \partial \alpha_2} \\ \frac { \partial u^3 } { \partial \alpha_3} \\ \end{array} \right) = L \cdot \left( \begin{array}{c} u_{10} \\ \frac { \partial u_{10} } { \partial \alpha_1} \\ u_{11} \\ \frac { \partial u_{11} } { \partial \alpha_1} \\ u_{12} \\ \frac { \partial u_{12} } { \partial \alpha_1} \\ u_{30} \\ \frac { \partial u_{30} } { \partial \alpha_1} \\ u_{31} \\ \frac { \partial u_{31} } { \partial \alpha_1} \\ u_{32} \\ \frac { \partial u_{32} } { \partial \alpha_1} \\ \end{array} \right) $


In [53]:
L=zeros(12,12)
h=Symbol('h')
p0=1/2-alpha3/h
p1=1/2+alpha3/h
p2=1-(2*alpha3/h)**2

L[0,0]=p0
L[0,2]=p1
L[0,4]=p2

L[1,1]=p0
L[1,3]=p1
L[1,5]=p2

L[3,0]=p0.diff(alpha3)
L[3,2]=p1.diff(alpha3)
L[3,4]=p2.diff(alpha3)

L[8,6]=p0
L[8,8]=p1
L[8,10]=p2

L[9,7]=p0
L[9,9]=p1
L[9,11]=p2

L[11,6]=p0.diff(alpha3)
L[11,8]=p1.diff(alpha3)
L[11,10]=p2.diff(alpha3)

L


Out[53]:
$$\left[\begin{array}{cccccccccccc}- \frac{\alpha_{3}}{h} + 0.5 & 0 & \frac{\alpha_{3}}{h} + 0.5 & 0 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & - \frac{\alpha_{3}}{h} + 0.5 & 0 & \frac{\alpha_{3}}{h} + 0.5 & 0 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{1}{h} & 0 & \frac{1}{h} & 0 & - \frac{8 \alpha_{3}}{h^{2}} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - \frac{\alpha_{3}}{h} + 0.5 & 0 & \frac{\alpha_{3}}{h} + 0.5 & 0 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & - \frac{\alpha_{3}}{h} + 0.5 & 0 & \frac{\alpha_{3}}{h} + 0.5 & 0 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - \frac{1}{h} & 0 & \frac{1}{h} & 0 & - \frac{8 \alpha_{3}}{h^{2}} & 0\end{array}\right]$$

In [64]:
B_General = zeros(9, 12)
B_General[0,1] = S(1)
B_General[1,2] = S(1)
B_General[2,3] = S(1)

B_General[3,5] = S(1)
B_General[4,6] = S(1)
B_General[5,7] = S(1)

B_General[6,9] = S(1)
B_General[7,10] = S(1)
B_General[8,11] = S(1)

for row_index in range(9):
    i,j=row_index_to_i_j_grad(row_index)
    
    B_General[row_index, 0] = -Symbol("G_{{{}{}}}^1".format(i+1,j+1))
    B_General[row_index, 4] = -Symbol("G_{{{}{}}}^2".format(i+1,j+1))
    B_General[row_index, 8] = -Symbol("G_{{{}{}}}^3".format(i+1,j+1))

B_General


Out[64]:
$$\left[\begin{array}{cccccccccccc}- G_{11}^1 & 1 & 0 & 0 & - G_{11}^2 & 0 & 0 & 0 & - G_{11}^3 & 0 & 0 & 0\\- G_{12}^1 & 0 & 1 & 0 & - G_{12}^2 & 0 & 0 & 0 & - G_{12}^3 & 0 & 0 & 0\\- G_{13}^1 & 0 & 0 & 1 & - G_{13}^2 & 0 & 0 & 0 & - G_{13}^3 & 0 & 0 & 0\\- G_{21}^1 & 0 & 0 & 0 & - G_{21}^2 & 1 & 0 & 0 & - G_{21}^3 & 0 & 0 & 0\\- G_{22}^1 & 0 & 0 & 0 & - G_{22}^2 & 0 & 1 & 0 & - G_{22}^3 & 0 & 0 & 0\\- G_{23}^1 & 0 & 0 & 0 & - G_{23}^2 & 0 & 0 & 1 & - G_{23}^3 & 0 & 0 & 0\\- G_{31}^1 & 0 & 0 & 0 & - G_{31}^2 & 0 & 0 & 0 & - G_{31}^3 & 1 & 0 & 0\\- G_{32}^1 & 0 & 0 & 0 & - G_{32}^2 & 0 & 0 & 0 & - G_{32}^3 & 0 & 1 & 0\\- G_{33}^1 & 0 & 0 & 0 & - G_{33}^2 & 0 & 0 & 0 & - G_{33}^3 & 0 & 0 & 1\end{array}\right]$$

In [66]:
simplify(B_General*L)


Out[66]:
$$\left[\begin{array}{cccccccccccc}\frac{G_{11}^1}{h} \left(\alpha_{3} - 0.5 h\right) & - \frac{\alpha_{3}}{h} + 0.5 & - \frac{G_{11}^1}{h} \left(\alpha_{3} + 0.5 h\right) & \frac{\alpha_{3}}{h} + 0.5 & \frac{4 G_{11}^1}{h^{2}} \alpha_{3}^{2} - G_{11}^1 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1 & \frac{G_{11}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{11}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{11}^3}{h^{2}} \alpha_{3}^{2} - G_{11}^3 & 0\\\frac{G_{12}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{12}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{12}^1}{h^{2}} \alpha_{3}^{2} - G_{12}^1 & 0 & \frac{G_{12}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{12}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{12}^3}{h^{2}} \alpha_{3}^{2} - G_{12}^3 & 0\\\frac{1}{h} \left(G_{13}^1 \left(\alpha_{3} - 0.5 h\right) - 1\right) & 0 & \frac{1}{h} \left(- G_{13}^1 \left(\alpha_{3} + 0.5 h\right) + 1\right) & 0 & \frac{1}{h^{2}} \left(G_{13}^1 \left(4 \alpha_{3}^{2} - h^{2}\right) - 8 \alpha_{3}\right) & 0 & \frac{G_{13}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{13}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{13}^3}{h^{2}} \alpha_{3}^{2} - G_{13}^3 & 0\\\frac{G_{21}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{21}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{21}^1}{h^{2}} \alpha_{3}^{2} - G_{21}^1 & 0 & \frac{G_{21}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{21}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{21}^3}{h^{2}} \alpha_{3}^{2} - G_{21}^3 & 0\\\frac{G_{22}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{22}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{22}^1}{h^{2}} \alpha_{3}^{2} - G_{22}^1 & 0 & \frac{G_{22}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{22}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{22}^3}{h^{2}} \alpha_{3}^{2} - G_{22}^3 & 0\\\frac{G_{23}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{23}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{23}^1}{h^{2}} \alpha_{3}^{2} - G_{23}^1 & 0 & \frac{G_{23}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{23}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{23}^3}{h^{2}} \alpha_{3}^{2} - G_{23}^3 & 0\\\frac{G_{31}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{31}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{31}^1}{h^{2}} \alpha_{3}^{2} - G_{31}^1 & 0 & \frac{G_{31}^3}{h} \left(\alpha_{3} - 0.5 h\right) & - \frac{\alpha_{3}}{h} + 0.5 & - \frac{G_{31}^3}{h} \left(\alpha_{3} + 0.5 h\right) & \frac{\alpha_{3}}{h} + 0.5 & \frac{4 G_{31}^3}{h^{2}} \alpha_{3}^{2} - G_{31}^3 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1\\\frac{G_{32}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{32}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{32}^1}{h^{2}} \alpha_{3}^{2} - G_{32}^1 & 0 & \frac{G_{32}^3}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{32}^3}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{32}^3}{h^{2}} \alpha_{3}^{2} - G_{32}^3 & 0\\\frac{G_{33}^1}{h} \left(\alpha_{3} - 0.5 h\right) & 0 & - \frac{G_{33}^1}{h} \left(\alpha_{3} + 0.5 h\right) & 0 & \frac{4 G_{33}^1}{h^{2}} \alpha_{3}^{2} - G_{33}^1 & 0 & \frac{1}{h} \left(G_{33}^3 \left(\alpha_{3} - 0.5 h\right) - 1\right) & 0 & \frac{1}{h} \left(- G_{33}^3 \left(\alpha_{3} + 0.5 h\right) + 1\right) & 0 & \frac{1}{h^{2}} \left(G_{33}^3 \left(4 \alpha_{3}^{2} - h^{2}\right) - 8 \alpha_{3}\right) & 0\end{array}\right]$$

In [70]:
simplify(E*B*D*P*L)


Out[70]:
$$\left[\begin{array}{cccccccccccc}0 & - \frac{1}{R h} \left(R + \alpha_{3}\right) \left(\alpha_{3} - 0.5 h\right) & 0 & \frac{1}{R h} \left(R + \alpha_{3}\right) \left(\alpha_{3} + 0.5 h\right) & 0 & - \frac{1}{R h^{2}} \left(R + \alpha_{3}\right) \left(4 \alpha_{3}^{2} - h^{2}\right) & - \frac{1}{R^{2} h} \left(R + \alpha_{3}\right) \left(\alpha_{3} - 0.5 h\right) & 0 & \frac{1}{R^{2} h} \left(R + \alpha_{3}\right) \left(\alpha_{3} + 0.5 h\right) & 0 & - \frac{1}{R^{2} h^{2}} \left(R + \alpha_{3}\right) \left(4 \alpha_{3}^{2} - h^{2}\right) & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - \frac{1}{h} & 0 & \frac{1}{h} & 0 & - \frac{8 \alpha_{3}}{h^{2}} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{1}{h} - \frac{0.5}{R} & 0 & \frac{1}{h} - \frac{0.5}{R} & 0 & \frac{1}{R h^{2}} \left(4 \alpha_{3}^{2} - 8 \alpha_{3} \left(R + \alpha_{3}\right) - h^{2}\right) & 0 & 0 & - \frac{\alpha_{3}}{h} + 0.5 & 0 & \frac{\alpha_{3}}{h} + 0.5 & 0 & - \frac{4 \alpha_{3}^{2}}{h^{2}} + 1\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$

In [73]:
D_p_L = StrainL*L
K = Symbol('K')
D_p_L = D_p_L.subs(R, 1/K)
simplify(D_p_L)


Out[73]:
$$\left[\begin{array}{cccccccccccc}0 & \frac{- \alpha_{3} + 0.5 h}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{\alpha_{3} + 0.5 h}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{- 4 \alpha_{3}^{2} + h^{2}}{h^{2} \left(K \alpha_{3} + 1\right)} & - \frac{K \left(\alpha_{3} - 0.5 h\right)}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{K \left(\alpha_{3} + 0.5 h\right)}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{K \left(- 4 \alpha_{3}^{2} + h^{2}\right)}{h^{2} \left(K \alpha_{3} + 1\right)} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & - \frac{1}{h} & 0 & \frac{1}{h} & 0 & - \frac{8 \alpha_{3}}{h^{2}} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{0.5 K h + 1.0}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{- 0.5 K h + 1.0}{h \left(K \alpha_{3} + 1\right)} & 0 & - \frac{4 K \alpha_{3}^{2} + K h^{2} + 8 \alpha_{3}}{h^{2} \left(K \alpha_{3} + 1\right)} & 0 & 0 & \frac{- \alpha_{3} + 0.5 h}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{\alpha_{3} + 0.5 h}{h \left(K \alpha_{3} + 1\right)} & 0 & \frac{- 4 \alpha_{3}^{2} + h^{2}}{h^{2} \left(K \alpha_{3} + 1\right)}\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$

In [77]:
p02=integrate((0.5-alpha3/h)*(1-(2*alpha3/h)**2), (alpha3, -h/2, h/2))
expand(simplify(p02))


Out[77]:
$$0.333333333333333 h$$

Virtual work


In [48]:
Q=E*B*D*P
S=simplify(Q.T*C_isotropic_matrix_alpha*Q*(1+alpha3/R)**2)
S


Out[48]:
$$\left[\begin{array}{cccccccccccc}\frac{\mu}{R^{2}} & 0 & 0 & - \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & - \frac{\mu}{R} & 0 & 0\\0 & \lambda + 2 \mu & 0 & 0 & 0 & 0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right) & 0 & \frac{1}{R} \left(\lambda + 2 \mu\right) & 0 & 0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right)\\0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0 & \mu & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} \left(\lambda + 2 \mu\right) & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right)^{2}\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0\\0 & \frac{1}{R} \left(\lambda + 2 \mu\right) & 0 & 0 & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right) & 0 & \frac{1}{R^{2}} \left(\lambda + 2 \mu\right) & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right)\\- \frac{\mu}{R} & 0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & \mu & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0\\0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} \left(\lambda + 2 \mu\right)\end{array}\right]$$

Isotropic material physical coordinates


In [49]:
S = simplify(D_p.T*C_isotropic_matrix_alpha_p*D_p*(1+alpha3/R)**2)
S


Out[49]:
$$\left[\begin{array}{cccccccccccc}\frac{\mu}{R^{2}} & 0 & 0 & - \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & - \frac{\mu}{R} & 0 & 0\\0 & \lambda + 2 \mu & 0 & 0 & 0 & 0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right) & 0 & \frac{1}{R} \left(\lambda + 2 \mu\right) & 0 & 0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right)\\0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0 & \mu & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} \left(\lambda + 2 \mu\right) & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right)^{2}\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0\\0 & \frac{1}{R} \left(\lambda + 2 \mu\right) & 0 & 0 & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right) & 0 & \frac{1}{R^{2}} \left(\lambda + 2 \mu\right) & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right)\\- \frac{\mu}{R} & 0 & 0 & \frac{\mu}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & 0 & \mu & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0\\0 & \frac{\lambda}{R} \left(R + \alpha_{3}\right) & 0 & 0 & 0 & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right)^{2} & 0 & \frac{\lambda}{R^{2}} \left(R + \alpha_{3}\right) & 0 & 0 & \frac{1}{R^{2}} \left(R + \alpha_{3}\right)^{2} \left(\lambda + 2 \mu\right)\end{array}\right]$$

Isotropic material physical coordinates - Tymoshenko


In [52]:
W = simplify(D_p_T.T*C_isotropic_matrix_alpha_p*D_p_T*(1+K*alpha3)**2)
W


Out[52]:
$$\left[\begin{matrix}K^{2} \mu & 0 & - K \mu & 0 & 0 & - K \mu\\0 & \lambda + 2 \mu & 0 & \alpha_{3} \left(\lambda + 2 \mu\right) & K \left(\lambda + 2 \mu\right) & 0\\- K \mu & 0 & \mu & 0 & 0 & \mu\\0 & \alpha_{3} \left(\lambda + 2 \mu\right) & 0 & \alpha_{3}^{2} \left(\lambda + 2 \mu\right) & K \alpha_{3} \left(\lambda + 2 \mu\right) & 0\\0 & K \left(\lambda + 2 \mu\right) & 0 & K \alpha_{3} \left(\lambda + 2 \mu\right) & K^{2} \left(\lambda + 2 \mu\right) & 0\\- K \mu & 0 & \mu & 0 & 0 & \mu\end{matrix}\right]$$

In [58]:
h=Symbol('h')
E=Symbol('E')
v=Symbol('nu')
W_a3 = integrate(W, (alpha3, -h/2, h/2))
W_a3 = W_a3.subs(la, E*v/((1+v)*(1-2*v))).subs(mu, E/((1+v)*2))
W_a3 = simplify(W_a3)
W_a3


Out[58]:
$$\left[\begin{matrix}\frac{E K^{2} h}{2 \left(\nu + 1\right)} & 0 & - \frac{E K h}{2 \nu + 2} & 0 & 0 & - \frac{E K h}{2 \nu + 2}\\0 & \frac{E h \left(\nu - 1\right)}{2 \nu^{2} + \nu - 1} & 0 & 0 & \frac{E K h \left(\nu - 1\right)}{2 \nu^{2} + \nu - 1} & 0\\- \frac{E K h}{2 \nu + 2} & 0 & \frac{E h}{2 \left(\nu + 1\right)} & 0 & 0 & \frac{E h}{2 \left(\nu + 1\right)}\\0 & 0 & 0 & \frac{E h^{3} \left(\nu - 1\right)}{12 \left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & 0\\0 & \frac{E K h \left(\nu - 1\right)}{2 \nu^{2} + \nu - 1} & 0 & 0 & \frac{E K^{2} h \left(\nu - 1\right)}{2 \nu^{2} + \nu - 1} & 0\\- \frac{E K h}{2 \nu + 2} & 0 & \frac{E h}{2 \left(\nu + 1\right)} & 0 & 0 & \frac{E h}{2 \left(\nu + 1\right)}\end{matrix}\right]$$

In [66]:
A_M = zeros(3)
A_M[0,0] = E*h/(1-v**2)
A_M[1,1] = 5*E*h/(12*(1+v))
A_M[2,2] = E*h**3/(12*(1-v**2))
Q_M = zeros(3,6)
Q_M[0,1] = 1
Q_M[0,4] = K
Q_M[1,0] = -K
Q_M[1,2] = 1
Q_M[1,5] = 1
Q_M[2,3] = 1

W_M=Q_M.T*A_M*Q_M
W_M


Out[66]:
$$\left[\begin{matrix}\frac{5 E K^{2} h}{12 \nu + 12} & 0 & - \frac{5 E K h}{12 \nu + 12} & 0 & 0 & - \frac{5 E K h}{12 \nu + 12}\\0 & \frac{E h}{- \nu^{2} + 1} & 0 & 0 & \frac{E K h}{- \nu^{2} + 1} & 0\\- \frac{5 E K h}{12 \nu + 12} & 0 & \frac{5 E h}{12 \nu + 12} & 0 & 0 & \frac{5 E h}{12 \nu + 12}\\0 & 0 & 0 & \frac{E h^{3}}{- 12 \nu^{2} + 12} & 0 & 0\\0 & \frac{E K h}{- \nu^{2} + 1} & 0 & 0 & \frac{E K^{2} h}{- \nu^{2} + 1} & 0\\- \frac{5 E K h}{12 \nu + 12} & 0 & \frac{5 E h}{12 \nu + 12} & 0 & 0 & \frac{5 E h}{12 \nu + 12}\end{matrix}\right]$$

Isotropic material physical coordinates - Square


In [61]:
W = simplify(D_p_L.T*C_isotropic_matrix_alpha_p*D_p_L*(1+K*alpha3)**2)
W


Out[61]:
$$\left[\begin{array}{cccccccccccc}\frac{\mu}{h^{2}} \left(0.5 K h + 1\right)^{2} & 0 & 0.25 K^{2} \mu - \frac{1.0 \mu}{h^{2}} & 0 & \frac{\mu}{h^{3}} \left(- K \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & - \frac{\mu}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right)\\0 & \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right)^{2} \left(\lambda + 2 \mu\right) & 0 & - \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) & 0 & \frac{1}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) & \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & - \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{1}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) & 0\\0.25 K^{2} \mu - \frac{1.0 \mu}{h^{2}} & 0 & \frac{\mu}{h^{2}} \left(0.5 K h - 1\right)^{2} & 0 & \frac{\mu}{h^{3}} \left(- K \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) \left(- K \alpha_{3} + K \left(\alpha_{3} + 0.5 h\right) - 1\right) & 0 & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(- K \alpha_{3} + K \left(\alpha_{3} + 0.5 h\right) - 1\right) & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \alpha_{3} - K \left(\alpha_{3} + 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(- K \alpha_{3} + K \left(\alpha_{3} + 0.5 h\right) - 1\right)\\0 & - \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) & 0 & \frac{1}{h^{2}} \left(\alpha_{3} + 0.5 h\right)^{2} \left(\lambda + 2 \mu\right) & 0 & - \frac{1}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) & - \frac{1}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{1}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & - \frac{1}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) & 0\\\frac{\mu}{h^{3}} \left(- K \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{3}} \left(- K \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) \left(- K \alpha_{3} + K \left(\alpha_{3} + 0.5 h\right) - 1\right) & 0 & \frac{\mu}{h^{4}} \left(4 K \alpha_{3}^{2} + K h^{2} + 8 \alpha_{3}\right)^{2} & 0 & 0 & \frac{\mu}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(- K \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{\mu}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) - 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{\mu}{h^{4}} \left(16 K \alpha_{3}^{4} - K h^{4} + 32 \alpha_{3}^{3} - 8 \alpha_{3} h^{2}\right)\\0 & \frac{1}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) & 0 & - \frac{1}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) & 0 & \frac{1}{h^{4}} \left(4 \alpha_{3}^{2} - h^{2}\right)^{2} \left(\lambda + 2 \mu\right) & \frac{1}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & - \frac{1}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{1}{h^{4}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) & 0\\0 & \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & - \frac{1}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{1}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & \frac{1}{h^{2}} \left(K \left(\alpha_{3} - 0.5 h\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) + \left(K \alpha_{3} + 1\right) \left(K \lambda \left(\alpha_{3} - 0.5 h\right) + \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0 & - \frac{1}{h^{2}} \left(K \left(\alpha_{3} + 0.5 h\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) + \left(K \alpha_{3} + 1\right) \left(K \lambda \left(\alpha_{3} - 0.5 h\right) + \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0 & \frac{1}{h^{3}} \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(\alpha_{3} - 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right) \left(K \lambda \left(\alpha_{3} - 0.5 h\right) + \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0\\\frac{\mu}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(- K \alpha_{3} + K \left(\alpha_{3} + 0.5 h\right) - 1\right) & 0 & \frac{\mu}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(- K \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) & 0 & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} - 0.5 h\right)^{2} & 0 & - \frac{\alpha_{3}^{2} \mu}{h^{2}} + 0.25 \mu & 0 & \frac{\mu}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right)\\0 & - \frac{1}{h^{2}} \left(\alpha_{3} - 0.5 h\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{1}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & - \frac{1}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) & - \frac{1}{h^{2}} \left(K \left(\alpha_{3} - 0.5 h\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) + \left(K \alpha_{3} + 1\right) \left(K \lambda \left(\alpha_{3} + 0.5 h\right) + \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0 & \frac{1}{h^{2}} \left(K \left(\alpha_{3} + 0.5 h\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) + \left(K \alpha_{3} + 1\right) \left(K \lambda \left(\alpha_{3} + 0.5 h\right) + \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0 & - \frac{1}{h^{3}} \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(\alpha_{3} + 0.5 h\right) \left(\lambda + 2 \mu\right) + \lambda \left(K \alpha_{3} + 1\right)\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right) \left(K \lambda \left(\alpha_{3} + 0.5 h\right) + \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0\\- \frac{\mu}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} + 0.5 h\right) \left(K \alpha_{3} - K \left(\alpha_{3} + 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) - 8 \alpha_{3} \left(K \alpha_{3} + 1\right)\right) & 0 & 0 & - \frac{\alpha_{3}^{2} \mu}{h^{2}} + 0.25 \mu & 0 & \frac{\mu}{h^{2}} \left(\alpha_{3} + 0.5 h\right)^{2} & 0 & - \frac{\mu}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right)\\0 & \frac{1}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & - \frac{1}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) & 0 & \frac{1}{h^{4}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) & \frac{1}{h^{3}} \left(K \left(\alpha_{3} - 0.5 h\right) \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) + \left(K \alpha_{3} + 1\right) \left(K \lambda \left(4 \alpha_{3}^{2} - h^{2}\right) + 8 \alpha_{3} \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0 & \frac{1}{h^{3}} \left(- K \left(\alpha_{3} + 0.5 h\right) \left(- K \left(- 4 \alpha_{3}^{2} + h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) + \left(K \alpha_{3} + 1\right) \left(K \lambda \left(- 4 \alpha_{3}^{2} + h^{2}\right) - 8 \alpha_{3} \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0 & \frac{1}{h^{4}} \left(K \left(4 \alpha_{3}^{2} - h^{2}\right) \left(- K \left(- 4 \alpha_{3}^{2} + h^{2}\right) \left(\lambda + 2 \mu\right) + 8 \alpha_{3} \lambda \left(K \alpha_{3} + 1\right)\right) + 8 \alpha_{3} \left(K \alpha_{3} + 1\right) \left(- K \lambda \left(- 4 \alpha_{3}^{2} + h^{2}\right) + 8 \alpha_{3} \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\right)\right) & 0\\\frac{\mu}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(K \alpha_{3} - K \left(\alpha_{3} - 0.5 h\right) + 1\right) & 0 & \frac{\mu}{h^{3}} \left(4 \alpha_{3}^{2} - h^{2}\right) \left(- K \alpha_{3} + K \left(\alpha_{3} + 0.5 h\right) - 1\right) & 0 & \frac{\mu}{h^{4}} \left(16 K \alpha_{3}^{4} - K h^{4} + 32 \alpha_{3}^{3} - 8 \alpha_{3} h^{2}\right) & 0 & 0 & \frac{\mu}{h^{3}} \left(\alpha_{3} - 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right) & 0 & - \frac{\mu}{h^{3}} \left(\alpha_{3} + 0.5 h\right) \left(4 \alpha_{3}^{2} - h^{2}\right) & 0 & \frac{\mu}{h^{4}} \left(4 \alpha_{3}^{2} - h^{2}\right)^{2}\end{array}\right]$$

In [62]:
h=Symbol('h')
E=Symbol('E')
v=Symbol('nu')
W_a3 = integrate(W, (alpha3, -h/2, h/2))
W_a3 = W_a3.subs(la, E*v/((1+v)*(1-2*v))).subs(mu, E/((1+v)*2))
W_a3 = simplify(W_a3)
W_a3


Out[62]:
$$\left[\begin{array}{cccccccccccc}\frac{E \left(0.5 K h + 1\right)^{2}}{2 h \left(\nu + 1\right)} & 0 & \frac{E \left(0.5 K^{2} h^{2} - 2.0\right)}{4 h \left(\nu + 1\right)} & 0 & \frac{E K \left(0.666666666666667 K h + 1.33333333333333\right)}{2 \left(\nu + 1\right)} & 0 & 0 & - \frac{E \left(0.25 K h + 0.5\right)}{2 \nu + 2} & 0 & - \frac{E \left(0.25 K h + 0.5\right)}{2 \nu + 2} & 0 & - \frac{E \left(0.333333333333333 K h + 0.666666666666667\right)}{2 \nu + 2}\\0 & \frac{0.333333333333333 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{0.166666666666667 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{0.333333333333333 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & \frac{1.0 E \left(1.0 K h \nu - 1.33333333333333 K h + 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(1.0 K h \nu - 0.666666666666667 K h - 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(2.66666666666667 K h \nu - 1.33333333333333 K h - 2.66666666666667 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0\\\frac{E \left(0.5 K^{2} h^{2} - 2.0\right)}{4 h \left(\nu + 1\right)} & 0 & \frac{E \left(0.5 K h - 1\right)^{2}}{2 h \left(\nu + 1\right)} & 0 & \frac{E K \left(0.666666666666667 K h - 1.33333333333333\right)}{2 \left(\nu + 1\right)} & 0 & 0 & \frac{E \left(- 0.25 K h + 0.5\right)}{2 \left(\nu + 1\right)} & 0 & \frac{E \left(- 0.25 K h + 0.5\right)}{2 \left(\nu + 1\right)} & 0 & \frac{E \left(- 0.333333333333333 K h + 0.666666666666667\right)}{2 \left(\nu + 1\right)}\\0 & \frac{0.166666666666667 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{0.333333333333333 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{0.333333333333333 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & \frac{1.0 E \left(1.0 K h \nu - 0.666666666666667 K h + 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(1.0 K h \nu - 1.33333333333333 K h - 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(2.66666666666667 K h \nu - 1.33333333333333 K h + 2.66666666666667 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0\\\frac{E K \left(0.666666666666667 K h + 1.33333333333333\right)}{2 \left(\nu + 1\right)} & 0 & \frac{E K \left(0.666666666666667 K h - 1.33333333333333\right)}{2 \left(\nu + 1\right)} & 0 & \frac{2 E \left(7 K^{2} h^{2} + 20\right)}{15 h \left(\nu + 1\right)} & 0 & 0 & \frac{E \left(- 0.666666666666667 K h + 0.666666666666667\right)}{2 \left(\nu + 1\right)} & 0 & - \frac{0.333333333333333 E \left(K h + 1\right)}{\nu + 1} & 0 & - \frac{2 E K h}{5 \nu + 5}\\0 & \frac{0.333333333333333 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{0.333333333333333 E h \left(\nu - 1\right)}{\left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{8 E h \left(\nu - 1\right)}{15 \left(2 \nu^{2} + \nu - 1\right)} & \frac{1.0 E \left(0.666666666666667 K h \nu - 0.666666666666667 K h + 1.33333333333333 \nu\right)}{4.0 \nu^{2} + 2.0 \nu - 2.0} & 0 & \frac{1.0 E \left(0.666666666666667 K h \nu - 0.666666666666667 K h - 1.33333333333333 \nu\right)}{4.0 \nu^{2} + 2.0 \nu - 2.0} & 0 & \frac{4 E K h \left(3 \nu - 2\right)}{15 \left(\nu + 1\right) \left(2 \nu - 1\right)} & 0\\0 & \frac{1.0 E \left(1.0 K h \nu - 1.33333333333333 K h + 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(1.0 K h \nu - 0.666666666666667 K h + 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(0.666666666666667 K h \nu - 0.666666666666667 K h + 1.33333333333333 \nu\right)}{4.0 \nu^{2} + 2.0 \nu - 2.0} & \frac{1.0 E}{h \left(8.0 \nu^{2} + 4.0 \nu - 4.0\right)} \left(1.0 K^{2} h^{2} \nu - 1.66666666666667 K^{2} h^{2} + 4.0 K h \nu + 4.0 \nu - 4.0\right) & 0 & \frac{1.0 E \left(1.0 K^{2} h^{2} \nu - 0.333333333333333 K^{2} h^{2} - 4.0 \nu + 4.0\right)}{h \left(8.0 \nu^{2} + 4.0 \nu - 4.0\right)} & 0 & \frac{E K \left(- 2.66666666666667 \nu + \left(2 \nu - 1\right) \left(0.666666666666667 K h + 2.66666666666667\right)\right)}{2 \left(\nu + 1\right) \left(2 \nu - 1\right)} & 0\\- \frac{E \left(0.25 K h + 0.5\right)}{2 \nu + 2} & 0 & \frac{E \left(- 0.25 K h + 0.5\right)}{2 \left(\nu + 1\right)} & 0 & \frac{E \left(- 0.666666666666667 K h + 0.666666666666667\right)}{2 \left(\nu + 1\right)} & 0 & 0 & \frac{0.166666666666667 E h}{\nu + 1} & 0 & \frac{0.0833333333333333 E h}{\nu + 1} & 0 & \frac{0.166666666666667 E h}{\nu + 1}\\0 & \frac{1.0 E \left(1.0 K h \nu - 0.666666666666667 K h - 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(1.0 K h \nu - 1.33333333333333 K h - 2.0 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(0.666666666666667 K h \nu - 0.666666666666667 K h - 1.33333333333333 \nu\right)}{4.0 \nu^{2} + 2.0 \nu - 2.0} & \frac{1.0 E \left(1.0 K^{2} h^{2} \nu - 0.333333333333333 K^{2} h^{2} - 4.0 \nu + 4.0\right)}{h \left(8.0 \nu^{2} + 4.0 \nu - 4.0\right)} & 0 & \frac{1.0 E}{h \left(8.0 \nu^{2} + 4.0 \nu - 4.0\right)} \left(1.0 K^{2} h^{2} \nu - 1.66666666666667 K^{2} h^{2} - 4.0 K h \nu + 4.0 \nu - 4.0\right) & 0 & \frac{E K \left(2.66666666666667 \nu + \left(2 \nu - 1\right) \left(0.666666666666667 K h - 2.66666666666667\right)\right)}{2 \left(\nu + 1\right) \left(2 \nu - 1\right)} & 0\\- \frac{E \left(0.25 K h + 0.5\right)}{2 \nu + 2} & 0 & \frac{E \left(- 0.25 K h + 0.5\right)}{2 \left(\nu + 1\right)} & 0 & - \frac{0.333333333333333 E \left(K h + 1\right)}{\nu + 1} & 0 & 0 & \frac{0.0833333333333333 E h}{\nu + 1} & 0 & \frac{0.166666666666667 E h}{\nu + 1} & 0 & \frac{0.166666666666667 E h}{\nu + 1}\\0 & \frac{1.0 E \left(2.66666666666667 K h \nu - 1.33333333333333 K h - 2.66666666666667 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{1.0 E \left(2.66666666666667 K h \nu - 1.33333333333333 K h + 2.66666666666667 \nu\right)}{8.0 \nu^{2} + 4.0 \nu - 4.0} & 0 & \frac{4 E K h \left(3 \nu - 2\right)}{15 \left(\nu + 1\right) \left(2 \nu - 1\right)} & \frac{E K \left(- 2.66666666666667 \nu + \left(2 \nu - 1\right) \left(0.666666666666667 K h + 2.66666666666667\right)\right)}{2 \left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{E K \left(2.66666666666667 \nu + \left(2 \nu - 1\right) \left(0.666666666666667 K h - 2.66666666666667\right)\right)}{2 \left(\nu + 1\right) \left(2 \nu - 1\right)} & 0 & \frac{4 E \left(7 K^{2} h^{2} \nu - 5 K^{2} h^{2} + 20 \nu - 20\right)}{15 h \left(2 \nu^{2} + \nu - 1\right)} & 0\\- \frac{E \left(0.333333333333333 K h + 0.666666666666667\right)}{2 \nu + 2} & 0 & \frac{E \left(- 0.333333333333333 K h + 0.666666666666667\right)}{2 \left(\nu + 1\right)} & 0 & - \frac{2 E K h}{5 \nu + 5} & 0 & 0 & \frac{0.166666666666667 E h}{\nu + 1} & 0 & \frac{0.166666666666667 E h}{\nu + 1} & 0 & \frac{4 E h}{15 \left(\nu + 1\right)}\end{array}\right]$$

Mass matrix in physical coordinates


In [67]:
rho=Symbol('rho')
B_h=zeros(3,12)
B_h[0,0]=1
B_h[1,4]=1
B_h[1,8]=1
M=simplify(rho*P.T*B_h.T*G_con*B_h*P)
M


Out[67]:
$$\left[\begin{array}{cccccccccccc}\rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \rho & 0 & 0 & 0 & \rho & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \rho & 0 & 0 & 0 & \rho & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$

In [74]:
M_p = rho*B_h.T*B_h*(1+alpha3/R)**2
M_p


Out[74]:
$$\left[\begin{array}{cccccccccccc}\rho \left(1 + \frac{\alpha_{3}}{R}\right)^{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \rho \left(1 + \frac{\alpha_{3}}{R}\right)^{2} & 0 & 0 & 0 & \rho \left(1 + \frac{\alpha_{3}}{R}\right)^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \rho \left(1 + \frac{\alpha_{3}}{R}\right)^{2} & 0 & 0 & 0 & \rho \left(1 + \frac{\alpha_{3}}{R}\right)^{2} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$

In [76]:
mass_matrix_func = lambdify((rho, R, alpha3), M_p, "numpy")
mass_matrix_func(100,10,20)


Out[76]:
array([[900.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0., 900.,   0.,   0.,   0., 900.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0., 900.,   0.,   0.,   0., 900.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.],
       [  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,
          0.]])

In [77]:
stiffness_matrix_func = lambdify([R, mu, la, alpha3], S, "numpy")
stiffness_matrix_func(100, 200, 300, 400)


Out[77]:
array([[ 2.00e-02,  0.00e+00,  0.00e+00, -1.00e+01,  0.00e+00,  0.00e+00,
         0.00e+00,  0.00e+00,  0.00e+00, -2.00e+00,  0.00e+00,  0.00e+00],
       [ 0.00e+00,  7.00e+02,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         1.50e+03,  0.00e+00,  7.00e+00,  0.00e+00,  0.00e+00,  1.50e+03],
       [ 0.00e+00,  0.00e+00,  5.00e+03,  0.00e+00,  0.00e+00,  1.00e+03,
         0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00],
       [-1.00e+01,  0.00e+00,  0.00e+00,  5.00e+03,  0.00e+00,  0.00e+00,
         0.00e+00,  0.00e+00,  0.00e+00,  1.00e+03,  0.00e+00,  0.00e+00],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00],
       [ 0.00e+00,  0.00e+00,  1.00e+03,  0.00e+00,  0.00e+00,  2.00e+02,
         0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00],
       [ 0.00e+00,  1.50e+03,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         1.75e+04,  0.00e+00,  1.50e+01,  0.00e+00,  0.00e+00,  7.50e+03],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         0.00e+00,  5.00e+03,  0.00e+00,  0.00e+00,  5.00e+03,  0.00e+00],
       [ 0.00e+00,  7.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         1.50e+01,  0.00e+00,  7.00e-02,  0.00e+00,  0.00e+00,  1.50e+01],
       [-2.00e+00,  0.00e+00,  0.00e+00,  1.00e+03,  0.00e+00,  0.00e+00,
         0.00e+00,  0.00e+00,  0.00e+00,  2.00e+02,  0.00e+00,  0.00e+00],
       [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         0.00e+00,  5.00e+03,  0.00e+00,  0.00e+00,  5.00e+03,  0.00e+00],
       [ 0.00e+00,  1.50e+03,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,
         7.50e+03,  0.00e+00,  1.50e+01,  0.00e+00,  0.00e+00,  1.75e+04]])

In [78]:
import fem.geometry as g
import fem.model as m
import fem.material as mat
import fem.solver as s
import fem.mesh as me
import plot


def generate_layers(thickness, layers_count, material):
    layer_top = thickness / 2
    layer_thickness = thickness / layers_count
    layers = set()
    for i in range(layers_count):
        layer = m.Layer(layer_top - layer_thickness, layer_top, material, i)
        layers.add(layer)
        layer_top -= layer_thickness
    return layers


def solve(width, curvature, thickness):
    layers_count = 1
    layers = generate_layers(thickness, layers_count, mat.IsotropicMaterial.steel())
    mesh = me.Mesh.generate(width, layers, N, M, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
#    geometry = g.CorrugatedCylindricalPlate(width, curvature, corrugation_amplitude, corrugation_frequency)
    geometry = g.CylindricalPlate(width, curvature)
#    geometry = g.Geometry()
    model = m.Model(geometry, layers, m.Model.FIXED_BOTTOM_LEFT_RIGHT_POINTS)
    return s.solve(model, mesh, stiffness_matrix, mass_matrix)

def stiffness_matrix(material, geometry, x1, x2, x3):
    return stiffness_matrix_func(1/geometry.curvature, material.mu(), material.lam(), x3)

def mass_matrix(material, geometry, x1, x2, x3):
    return mass_matrix_func(material.rho, 1/geometry.curvature, x3)


# r=2
# width = r*2*3.14
# curvature = 1/r

width = 2
curvature = 0.8
thickness = 0.05

N = 100
M = 10


results = solve(width, curvature, thickness)
results_index = 0
plot.plot_init_and_deformed_geometry(results[results_index], 0, width, -thickness / 2, thickness / 2, 0)
#plot.plot_init_geometry(results[results_index].geometry, 0, width, -thickness / 2, thickness / 2, 0)
# plot.plot_strain(results[results_index], 0, width, -thickness / 2, thickness / 2, 0)


to_print = 20
if (len(results) < to_print):
    to_print = len(results)

for i in range(to_print):
    print(results[i].freq)


/anaconda3/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
25.518397195470808
37.65055553043148
53.02521516884746
59.12459462472745
73.89906602067028
76.7238220818741
92.86152431436375
95.74961566522097
111.69858911730539
115.37879991578257
128.43308532025998
133.39281956188958
147.98227182643734
148.02904662879308
161.91555797565047
166.73482225476707
176.81831442359024
181.56173242142907
193.8669205843182
194.17303269234304

In [ ]: